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Abstract.  A  new  approach  to  DSMC  collision  modelling,  called  viscosity-DSMC  or  /i-DSMC,  is  described  in  which  the 
time-averaged  temperature  is  used  to  set  the  characteristic  collision  cross-section  in  each  cell  such  that  the  Chapman-Enskog 
viscosity  is  that  given  by  any  desired  viscosity  law  p  =  p  (T),  including  a  curve  fit  to  experimental  data.  For  example,  a  hard 
sphere  collision  model,  with  hard  sphere  collision  probability,  used  with  a  different  molecular  size  in  each  cell  can  reproduce 
a  Sutherland  viscosity  law.  Similarly,  a  variable  hard  sphere  collision  model  can  reproduce  the  viscosity  given  by  the  more 
complicated  generalized  hard  collision  model,  by  making  the  reference  cross-section  a  function  of  the  temperature.  This 
model  is  used  to  calculate  the  structure  of  a  plane  ID  shock  and  the  results  agree  closely  with  those  from  standard  DSMC 
using  the  GHS  model.  A  particularly  simple  method  is  to  use  the  Maxwell  VHS  model,  in  which  all  collision  pairs  are  equally 
likely,  to  produce  any  desired  viscosity  law.  The  time-averaged  cell  temperature  is  available  in  standard  DSMC  as  part  of  the 
procedures  which  determine  the  steady  state  flow  and  the  new  methods  are  as  fast  as,  or  faster  than  standard  DSMC.  Unlike 
more  complicated  models  with  realistic  viscosities,  the  new  procedures  are  compatible  with  the  Borgnakke-Larsen  energy 
exchange  scheme  and  the  established  chemistry  models  for  DSMC. 


VISCOSITY  OF  REAL  GASES  AND  THE  VHS  MODEL 

The  most  commonly  used  DSMC  collision  model  is  the  variable  hard  sphere  (VHS)  which  gives  rises  to  a  power 
law  viscosity  p  oc  T®.  Although  this  viscosity  law  is  reasonably  accurate  over  a  limited  range  of  temperature  for 
a  particular  gas  it  does  not  represent  a  realistic  viscosity  over  all  temperatures.  For  T  <  2000  K,  the  viscosity  of  a 
typical  gas  displays  a  significant  deviation  from  a  simple  power  law  as  a  result  of  the  long-range  attractive  forces 
between  molecules.  Fig.  1  shows  the  viscosity  of  argon  compared  with  the  two  extremes  of  the  power  law,  co  = 
corresponding  to  hard  spheres  with  constant  cross-section,  and  co  =  1  corresponding  to  a  ‘Maxwell  molecule’.  Also 


T/Ti 

FIGURE  1.  Viscosity  for:  — ,  Sutherland,  Eq.  1 ;  — ,  hard  sphere  p / p\  =  (T /T\ ) 5 ; - ,  linear  law  p/p\  =  T /T\ .  o  argon  data 

Kestin  etal.  [1]  T\  =TS  =  142  K,  px  =  1.129  x  1(T5  Nm"1  s"1. 
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shown  is  the  Sutherland  viscosity, 


H/m  =  (T/T^2  (l+TjT^/il+TjT),  (1) 

where  7^  =  142  K  and  jli=  jl(T\)  which  fits  the  data  better  than  any  power  law. 

When  using  the  VHS  model  in  any  particular  DSMC  calculation  one  can  determine  the  range  of  temperatures 
expected  in  the  flow  and  choose  VHS  parameters  which  match  the  viscosity  reasonably  well.  However  it  would  be 
better  to  have  a  collision  model  which  fitted  the  experimental  data  over  all  temperatures.  Realistic  potentials,  such  as 
the  Morse  potential  [2],  the  Lennard-Jones  potential  [3,  4]  and  the  Maitland-Smith  potential  [5]  have  been  used  in 
DSMC  as  have  other  collision  models  [6,  7]  which  produce  realistic  viscosity  laws,  but  these  are  not  widely  used.  This 
is  partly  because  of  their  relative  complexity,  but  also  because  of  the  difficult  of  using  the  Borgnakke-Larsen  energy 
exchange  scheme  [8]  when  the  collision  probability  does  not  match  that  of  the  VHS  collision  model.  The  combination 
of  the  Borgnakke-Larsen  (BL)  exchange  model  and  the  VHS  model  is  now  the  de  facto  standard  of  DSMC.  Similarly, 
the  DSMC  procedures  for  chemically  reacting  flow  [9]  are  built  around  the  VHS  model. 

Here  I  show  how  an  arbitrary  viscosity  law  fi  =  li(T)  can  be  realized  in  DSMC  using  simple  collision  models 
which  are  compatible  with  the  BL  exchange  scheme.  This  is  good  not  only  in  itself,  but  also  because  it  makes  the 
construction  of  hybrid  DSMC/Navier-Stokes  solvers  easier;  the  viscosity  law  in  both  the  DSMC  and  Navier-Stokes 
code  can  be  the  best  available,  rather  than  that  dictated  by  computational  practicality.  The  new  method  runs  as  fast,  or 
faster  than  the  standard  VHS  collision  model. 

The  new  method,  called  /I- DSMC,  is  to  adjust  the  size  of  any  simple  collision  model  based  on  the  local  temperature 
to  produce  any  desired  viscosity  at  that  temperature.  We  use  the  hard  sphere  and  VHS  collision  models  as  the  basis 
of  the  new  method  and  in  each  case  produce  a  viscosity  law  different  from  the  usual  one  for  that  collision  model. 
The  method  has  been  tested  in  Couette  flow  and  for  highly  non-equilibrium  flow  in  the  interior  of  a  shock  and  has 
been  shown  to  produce  essentially  the  same  results  as  DSMC  using  more  complicated  collision  models.  A  particular 
form  of  the  general  method,  based  on  the  Maxwell  limit  of  VHS,  is  described  and  compared  with  standard  DSMC  in 
a  zero-dimensional  velocity  relaxation  calculation  and  the  supersonic  flow  around  a  blunt-faced  cylinder.  In  all  cases 
the  new  method  agrees  well  with  standard  DSMC;  the  deviations  are  small  and  confined  to  small  regions  of  the  flow. 
In  all  cases  shown  here  a  monatomic  gas,  with  y  =  5/3  has  been  used. 


HARD  SPHERE  MODEL  WITH  SUTHERLAND  VISCOSITY 


The  Chapman-Enskog  viscosity  for  hard  spheres,  with  diameter  d  is 

_  5m  (7 zRT)l/2 
M-T6  a  ’ 

where  a  =  7id2  is  the  the  total  collision  cross-section  and  R  =  k/m  is  the  ordinary  gas  constant.  A  given  viscosity  law 
\i  =  \l(T)  can  be  implemented  by  setting  the  cross-section  in  each  cell  as 


5m  (7 rtf  f) 1/2 


16  n(T) 


(2) 


where  T  is  the  time-averaged  translational  kinetic  temperature  in  the  cell. 

This  method  has  been  tested  for  supersonic  Couette  flow  in  Ref.  [10],  using  the  Sutherland  viscosity  (Eq.  1)  for 
ll  in  Eq.  2.  The  measured  shear  stress  Tmeas  =  pcxcy ,  was  determined  from  the  steady-state  velocity  distribution  and 
compared  with  the  theoretical  shear  stress  Tt  =  fisuthdux/dy,  where  /lsuth  is  the  Sutherland  viscosity  evaluated  for  the 
cell  temperature.  The  ratio  Tmeas/Tt  in  Fig.  2,  is  close  to  unity  across  most  of  the  flow  except  in  the  Knudsen  layer 
near  the  wall.  The  figure  also  shows  that  /lmeas/ Ak  is  significantly  different  from  unity  when  /it  is  calculated  using  the 
hard  sphere  viscosity;  in  other  words  the  hard  sphere  and  Sutherland  viscosities  are  significantly  different  for  these 
temperatures.  Thus,  the  hard  sphere  collision  model  with  a  different  collision  cross-section  in  each  cell  displays  the 
Sutherland  law  as  required. 
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FIGURE  2.  Measured  shear  stress  in  Couette  flow  simulation  using  /i-DSMC  and  Sutherland  viscosity,  compared  with  theoretical 
shear  stress  Tt  =  (TmQas)dux/dymQas.  o,  /rt  from  Eq.  1  (Sutherland);  +,  pt  =  ( T /T\)1^2  (hard  sphere).  Origin  of  y-axis  mid-way 
between  the  plates.  Distance  between  plates  2 H.  Nominal  Knudsen  number  2p\/ (p\c\H)  =  0.01,  c  =  (%RT\/ 7i)1/2 ,  subscript  1 
denotes  undisturbed  gas  state.  Wall  temperature  TW  =  T\.  Wall  speed  Vw  =  3  ( 2RT\ )^2.  Only  1  in  5  points  shown.  From  Ref.  [11]. 


VHS  COLLISION  MODEL  WITH  GHS  VISCOSITY 

The  VHS  collision  model  has  a  total  collision  cross-section  given  by 

<y  =  Or{gr/g?V 

where  g  is  the  collision  speed,  or  is  a  reference  cross-section,  gr  is  a  constant  reference  speed  and  V  is  a  constant  in 
the  range  0  to  1/2.  With  isotropic  scattering  and  this  total  cross-section  the  theoretical  viscosity  is 

15m  (7lRT)l/2(4RT)v  _  15 m  (: kRT)1/2 

M  “  8r(4-u)  <7rg2t>  “  8r(4-u)  Or 

for  gr  =  (4 RT)1/2.  To  achieve  an  arbitrary  viscosity  law  fi  =  n{T)  the  reference  cross-section  in  each  cell  is  set  as 

...  15m  (kRT)1/2 

Gr  ^  “  8r(4-u)  II  {T) 

In  Ref.  [10]  this  method,  with  0  =  1/6  was  used  to  calculate  the  structure  of  a  plane  normal  shock,  with  the  viscosity 
specified  as  that  for  the  generalized  hard  sphere  (GHS)  [6] 

157T1/2  (T/7b)1/2+Ul  mgp 

M  16r(4-ui)  [0  +  (1-0)5]  (70  ‘ 

With  (j>  =  0.61,  0|  =2/13,  v2  =  14/13,  a0  =  6.457  x  10-19m2,  T0  =  300  K,  %  S0 (7b/r)U2_Ul  and  =  (47?7b)1/2, 
Eq.  3  hts  the  data  for  argon  better  than  the  Sutherland  law.  The  value  of  So  was  set  so  that  jd  (300)  =  2.272  x  10  Nm-1. 
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FIGURE  3.  (a)  Normalized  density  and  temperature  profiles  for  normal  shock  ( M\  =  4,  T\  =  150  K,  y  =  5/3).  Solid  lines  show 

DSMC  (GHS)  results;  symbols  are  for  modified  fi- DSMC.  o,  Tx ;  V,  Tp ;  o,  p.  (b)  Differences  between  Tx  calculated  by  DSMC 
(GHS)  and  p  -DSMC  for  the  shock  in  part  (a).  The  ‘error’  is  expressed  as  a  percentage  of  the  temperature  jump  across  the  shock,  o, 
standard  p-DSMC,  maximum  error  4.6%;  +,  modified  p-DSMC,  maximum  error  2.6%.  From  Ref.  [10] 


The  /i -DSMC  results  for  a  Mach  4  shock  with  this  viscosity  are  compared  in  Fig.  3  with  DSMC  results  using  the 
GHS  model1  for  which  the  total  collision  cross-section  is  given  by 

ct/cto  -  <t>  {go/gfVl  +  (1-<P)  (go/gf^2  ■ 

Fig.  3(a)  shows  the  shock  profiles  of  p,  Tx  and  parallel  temperature  Tp  =  (Ty  +  TZ)/ 2  for  both  methods.  The  new 
method  produces  results  almost  identical  to  those  for  standard  DSMC  using  the  GHS  collision  model.  There  is  a  slight 
difference  in  the  Tx  profiles  in  a  small  region  upstream  of  the  shock.  In  the  case  shown  a  small  modification  was  used 
for  p-DSMC;  the  size  of  collision  cross-section  was  calculated  from  Eq.  3  using  Tx ,  rather  than  T .  Fig.  3(b)  shows  the 
differences  in  the  Tx  profiles,  between  DSMC  and  p-DSMC,  the  latter  using  both  T  and  Tx .  The  maximum  deviation 
from  DSMC  is  4.6%  of  the  temperature  rise  across  the  shock  for  standard  p  -DSMC  and  2.6%  for  the  modified  method. 


MAXWELL  CROSS-SECTION  WITH  VISCOSITY  ~  7® 

The  ‘Maxwell  VHS’  model  [13]  is  the  special  case  (v  =  1/2)  of  the  VHS  model  for  which 

<7  =  Orgr/g. 

The  collision  probability  (<*  go)  is  independent  of  relative  speed  g  for  this  model.  The  collision  loops  are  particularly 
simple  which  can  increase  the  computational  speed.  When  the  reference  cross-section  in  each  cell  is  adjusted  according 
to  the  desired  viscosity  law,  the  collision  rate  and  number  of  collisions  required  in  one  time  step  At  become 

v  (T )  =  2  hkT  and  Nco\\  =  hkT  /fix  NAt 

respectively,  where  n  is  the  time-averaged  number  density  and  N  is  the  number  of  simulator  particles  in  the  cell.  Note 
that  the  use  of  n  is  standard  in  DSMC  [14]  and  that  T /p  (f )  can  be  calculated  at  regular  intervals  for  each  cell. 

This  method  was  called  ‘collision  rate  DSMC5  (v-DSMC)  in  Ref.  [11]  where  its  behavior  in  a  zero-dimensional 


1  A  slight  modification  of  the  collision  cross-section  for  g  <  go  was  used  which  makes  the  GHS  model  computationally  efficient,  with  negligible 
effect  on  its  theoretical  viscosity.  See  Refs.  [12]  and  [10]. 
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FIGURE  4.  Thermal  speed  distribution  after  one  nominal  collision  time  t  —  2p/  (pc2),  where  c  =  (8 RTe/7l)1^2,  is  the  mean 
thermal  speed  and  Te  is  the  equilibrium  (or  kinetic)  temperature.  The  initial  distribution  was  bi-modal  corresponding  to  a  high 
speed,  high  temperature  jet  mixed  with  a  low  temperature  gas  at  rest.  — ,  DSMC  (VHS)  with  v  =  0.31;  +,  v-DSMC  with  matching 
viscosity  p  (Te).  The  collision  rate  for  v-DSMC  «  21%  greater  than  for  DSMC  (VHS).  From  Ref.  [11] 


relaxation  calculation  was  compared  with  standard  DSMC.  Every  DSMC  calculation  requires  such  a  relaxation 
calculation  in  each  cell  at  each  time  step.  If,  starting  from  any  non-equilibrium  distribution  of  velocities,  two  different 
collision  models  give  the  same  (statistical)  velocity  distribution  at  the  end  of  the  collision  phase  then  the  two  models 
will  produce  identical  results  when  used  in  DSMC,  regardless  of  how  many  collisions  each  model  requires. 

In  Ref.  [11],  V-DSMC  was  tested  against  DSMC  for  a  number  of  highly  non-equilibrium  initial  velocity  distribu¬ 
tions.  An  example  is  shown  in  Fig.  4.  In  the  initial  state  1/5  of  the  particles  had  velocities  selected  from  an  equilibrium 
distribution  with  a  mean  (20,0,0)  and  a  characteristic  thermal  speed  of  3  (arbitrary  units).  The  remaining  4/5  had  a 
mean  velocity  (0,0,0)  and  characteristic  speed  of  1 .  The  thermal  speed  distributions  for  the  two  methods  are  virtually 
the  same  after  the  same  elapsed  time  t  =  2p/  (pc2) ,  where  c  =  (8 RT / n )^2  is  the  mean  thermal  speed.  The  collision 
rate  for  v-DSMC  was  about  21%  greater  than  for  standard  VHS;  a  larger  number  of  collision  is  required  in  v-DSMC 
because,  with  the  standard  VHS  model,  collision  pairs  are  weighted  towards  higher  relative  velocities  so  that  an  aver¬ 
age  VHS  collision  is  more  ‘efficient’  in  redistributing  energy.  This  same  effect  can  be  shown  for  the  Krook  and  Wu 
near-equilibrium  exact  solutions  of  the  Boltzmann  equation  [15];  different  differential  collision  cross-sections  give 
identical  results  if  their  (integrated)  viscosity  cross-sections  are  matched. 

Simulations  using  v-DSMC  and  VHS,  with  the  same  viscosity  law  p  T° 12  (v  =  0.22),  were  performed  for 
the  Mach  10  flow  around  a  blunt-ended  cylinder  with  its  blunt  face  normal  to  the  freestream  velocity.  The  nominal 
Knudsen  number  was  Kn^o o  =  2 p^/  (pooCooD)  =  0.03,  and  the  wall  temperature  ratio  Tw/Too  was  0.26.  Fig.  5(a)  shows 
Mach  number  contours  which  are  virtually  indistinguishable,  except  for  the  slight  difference  ahead  of  the  shock.  The 
agreement  in  the  expansion  round  the  sharp  corner  is  very  close.  Fig.  5(b)  shows  the  near  stagnation  line  profiles  of 
density  p  and  x-component  of  translation  temperature  Tx.  The  density  rises  through  the  shock  and  near  the  cold  wall 
where  the  temperature  falls.  The  steep  rise  of  Tx  in  the  shock,  ahead  of  the  density  rise,  can  be  seen.  The  v-DSMC 
and  DSMC  results  agree  except  for  the  small  difference  in  Tx  ahead  of  the  shock  which  was  also  seen  in  the  ID  shock 
profiles  of  Fig.  3. 


DISCUSSION  AND  CONCLUSIONS 

Although  it  has  not  been  demonstrated  here,  it  is  clear  that,  since  a  standard  VHS  collision  model  is  used  in  each  cell, 
the  new  method  is  compatible  with  the  Borgnakke-Farsen  energy  exchange  scheme  and  with  the  standard  chemistry 
models  of  DSMC,  which  are  also  based  on  the  VHS  model.  The  new  method  makes  the  VHS  model  display  an  arbitrary 
viscosity  relationship  by  adjusting  the  number  of  collisions  calculated  in  each  cell,  based  on  the  cell  temperature.  There 
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FIGURE  5.  Antisymmetric  flow  about  a  blunt-faced  cylinder.  v-DSMC  compared  with  standard  VHS  model.  Moo  =  10,  Tw /Too  = 
0.26,  r  =  5/3,  n/n oc  =  (T/T^fJ2.  KnD  ca  =  2p^/ (p^c^D)  =  0.03.  (a)  Mach  number  contours:  -  -,  VHS;  — ,  v-DSMC.  (b) 
Temperature  Tx  and  density  p  in  cells  adjacent  to  the  stagnation  streamline.  □,  VHS;  •,  v-DSMC. 


is  virtually  no  extra  computational  overhead  involved  in  obtaining  the  time-averaged  temperature  which  is  available 
in  DSMC  codes  as  the  evolution  of  the  flow  is  tracked  to  steady  state.  Although  the  cell  temperature  is  not  routinely 
used  in  DSMC  as  part  of  the  simulation  procedure,  Bird  has  suggested  [14]  that  it  could  be  used  to  vary  the  exchange 
factor  in  the  BL  exchange  scheme  to  agree  with  experimental  results. 

Any  value  of  v  in  the  VHS  model  can  be  used  as  the  basis  of  /r- DSMC.  For  v  =  0  (hard  sphere)  the  number  of 
collisions  calculated  at  each  step  is  as  low  as  possible,  while  for  v  =  1/2  (Maxwell  VHS)  the  collision-pair  selection 
is  the  easiest  and  the  new  method  can  be  twice  as  fast  as  DSMC,  although  this  depends  on  details  of  the  code  and  flow. 
For  the  sake  of  a  more  realistic  distribution  of  relative  velocities  in  collisions  an  intermediate  value  of  v  =  1  / 4  might 
be  better  than  either  of  the  extreme  values. 

The  new  method  has  been  shown  to  be  accurate  for  the  highly  non-equilibrium  flow  in  the  interior  of  a  shock  and 
for  the  high  speed  flow  around  a  blunt  cylinder  which  contains  an  expansion  around  the  sharp  corner  of  the  front  face. 
Although  more  testing  is  desirable,  the  new  method  appears  to  be  very  promising. 
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